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ABSTRACT 
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■ We present a new Eulerian code able to follow the evolution of large-scale structures in the 

Universe in the weakly nonlinear regime. We compare test results with a N-body code and analytical 
results. 
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1. Introduction 

The Universe is homogeneous on scales larger than at least 100 Mpc. There- 
fore, any simulations intended to completely predict the development of structures 
ought to work with a dynamical range 10 kpc-100 Mpc. A disadvantage of this fact 
is that enormous catalogues of galaxies are needed for a comparison of models of 
large-scale structure versus observations. The first sufficient catalogue will be the 
SDSS (Gunn and Knapp 1993), to come within a few years. On the other hand, a 
great advantage of examining very large scale structures is that density contrasts on 
these scales are relatively low (i.e., weakly nonlinear, with 5 = < 1), and var- 
ious fields observed (convolved) with large windows can be worked on either with 
the use of analytical methods (like perturbation theory) or numerical simulations 
with simplified physics (eg., effects of pressure, radiation etc can be neglected). 
Even in the weakly nonlinear regime there are many features which have not been 
examined with analytical methods, e.g., moments of projected velocities or redshift 
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distortions. All of the above give a good reason for writing a new, fast code, able 
to follow the evolution of structure in various cosmological models. 

We present a new Eulerian code, which can be applied to linear and weakly 
non-linear evolutionary stages of dense regions of the Universe and to trace the 
voids (iV-body codes are not good for describing low density regions). The code 
is tested and the results are compared with theoretical predictions, as well as with 
results obtained from modified AP 3 M N -body code by Couchman (1991), (1994). 



2. Numerical Method 



2.1. Description 

Cosmological Pressureless Parabolic Advection (CPPA) is an Eulerian code, 
solving the dynamical equations of a self-gravitating pressureless fluid (Poisson, 
continuity and Euler). 

<9 2 $ 

= A-kGq (1) 



dr l dr 
dQ + dgv^ = Q 
dr dr k 
dgu 1 dgu l u k <9<& 
dr dr k ® dr % 

The code works with a fixed and uniform 3-D Cartesian grid in cosmological 
frame. Periodic boundary conditions are employed, which do not introduce any 
systematic deviations at grid boundaries, and are easy to implement for solving 
Poisson equation using Fourier methods. The Poisson solver is adapted from a 
standard Helmholtz equation solving routine available in FISHPAK library (Adams, 
Schwarztrauber, and Sweet 1980). 

The advection algorithm is based on the PPM method of Colella and Woodward 
(1984) with pressure terms set to zero. We replaced the nonlinear Riemann solver 
by a simple procedure which defines the interface values as the upwind value equal 
to one of the effective states averaged over domains of dependence in neighboring 
zones. Since there is no sufficient justification for steepening of the density jumps 
we did not use PPM steepening algorithm. Also, the whole set of dissipation pro- 
cedures is not needed for present application. The temporal evolution is accurate 
to second order while the spatial advection is formally third order accurate. The 
non-cosmological version of the code has been widely tested and proved isotropic 
and precisely solving Euler and Poisson equations. 

To adapt the code for cosmological purposes we changed the variables and 
coordinates in Eq. (l)-(3) following the recipe given by Gnedin (1995). Let r, 
r l , u l , g, <I> denote time, positions, velocities, density and gravitational potential 
in physical units. We define new set of variables t, x l , v l , p and 4> given by 
Eqs (4)-(8): 
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where a is the scale factor, Ho is the Hubble constant today, VLq is the dimension- 
less density parameter today, G denotes the gravitational constant and L is a di- 
mensional parameter denoting the computational box size. The density p has been 
scaled to (p) = 1, so p = 1 + S. In these coordinates, assuming that scale factor 
changes with time according to Friedmann-Robertson- Walker solution, equations 
Eq. (l)-(3) are 



d 2 <j) 3aQ 



dx l dx l 2 



dp d(pv 



5, (9) 



and 

8{pjf) | djptfv 11 ) | p d(j) _ Q 
dt dx k dx % 

The form of Eqs (9)-(ll) is very similar to Eqs (l)-(3), expressed in physical coor- 
dinates, which makes it relatively easy to adapt any hydrodynamical code to work 
in the cosmological frame. For a given cosmological model, only the functions 
a(t) and r(t) have to be evaluated numerically before the simulation is started. 
The present version works with A = , but it can be easily modified to deal with 
models with a non-zero cosmological constant, which requires only adding one 
more term to Eq. (9). 

2.2. Initial Conditions 

For cosmological simulations we assume a Gaussian distribution of the initial 
density contrast, 5 , with a given power spectrum and amplitude in the linear regime 
((5 2 )5 < 0.05). First, the initial density field S(x) is set up in the form of its 
Fourier transform S(k) , with amplitudes according to a given power spectrum and 
phases chosen randomly. To keep the density field real we impose 

5(-k) = 5(h)* . (12) 

To prevent exceeding the Nyquist frequency, the spectra are multiplied by e ~( fe / fccut ) 
where k cn t is smaller than kj^q (typically, we used k cut /kN q = 0.5 . . . 0.7). Next, 
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the field is Fourier-transformed back to the real space. Simulations of cosmologi- 
cal large-scale structure are started at early stages, when the linear approximation 
in the perturbation theory (see e.g., Peebles 1980, §11) can be applied. According 
to this approximation 

S = A 1 D 1 {t)+A 2 D 2 (t), (13) 

where D\ and D 2 denote the growing and decaying modes respectively. We as- 
sume that our starting evolutionary stage is advanced enough for not taking the 
decaying mode into account, i.e., ^(fstart) <C D i(r s tart)- From the continuity 
Eq. (2) we find that the initial velocity u and density contrast 5 are then related as 
follows: 

^ + lv-n = 0. (14) 
dr a 

Hence, 

V-u = -aH/<5, (15) 

where H is the current value of the Hubble's expansion parameter and / is given 
by 

a d£>i _ a dDi 1 1 Di 
* ~ lTi~d^ ~ D^~dr da/dr ~ 
and can be well approximated by 

/K!i o, + ^( 1+ l n ). 

Q denotes here the density parameter at r star t , and S7a = A/3H 2 is related to the 
cosmological constant. The vector field u is then the gradient of a scalar field ip , 
defined by the Poisson equation 

Ay> = aH/<5. (16) 

Therefore, given the density, the Poisson solver can be used for evaluating the initial 
velocity field. 

3. TheAP 3 MCode 

Our comparison code is the adaptive P 3 M developed by Couchman (1991),(1994). 
In this paper we will note only our changes to the code. Apart of altering the pro- 
gram's I/O we made two major changes. First of them was to introduce a CFL type 
evaluation of next time-step length. The second one was to change the time coor- 
dinate as explained below, and adapt the code to work in cosmological frame with 
Q / 1 , which was not supported by the original version. In comoving coordinates, 
the equation of motion of i -th particle in the expanding Universe reads 
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where r denotes the physical time and a the scale factor. As in CPPA, we will use 
another time variable, t, defined by Eq. (4). Eq. (17) will transform to: 

d 2 Xj _ 3f? L 3 <sr^ Xj- xj 

dt 2 " a 8itN ^ \xi-Xj\z' { ' 

where L is the computational box size and all of the N particles have mass equal 

8 ° G ° N , which is taken to be the mass unit. It is worth mentioning that by intro- 
ducing t instead of physical time r we have got rid of the cosmological drag term 
(in the original Couchman's version of the code time-stepping was performed at 
constant-length steps in a a , where a is an arbitrary real exponent, and the calcu- 
lations were time-staggered due to presence of the drag term). 



4. Code Tests 



The CPPA code has been tested versus the modified AP 3 M and analytical re- 
sults (the latter including both exact solutions and approximate calculations of the 
perturbation theory). 

4.1. Zel'dovich Pancake Test 

Zel'dovich Pancake is the final stage of the evolution of a single-wave perturba- 
tion, resulting in formation of a caustic. It may be also regarded as an idealization 
of a wall between two voids. The chief advantage of the pancake problem is that 
for Oo = 1 there exists an exact analytical solution: 

q(x, a) = R ( 1 — — cos(fcx) 
V a c 

where a c is a parameter (the value of a at which an infinitely dense caustic is 

I 2 2 

formed), and R = ^ ac n - is the normalization factor. We performed this test on a 
64 3 grid, setting up 1, 2, 4, 8 and 16 waves with initial dispersion 0.003 at an epoch 
corresponding to the scale factor equal to 0.001 of its present value. 1-D sections 
of the initial conditions are presented in Fig. 1. Fig. 2 shows the resulting pancakes 
at a = 0.1. 

The density distribution in the AP 3 M code is represented by the particle coor- 
dinates, actually it is a sum of N Dirac delta functions centered on each of the mass 
tracers. To obtain the grid-based density field, we convolved it with a Triangular 
Shaped Cloud (TSC) kernel function 

U{x,y,z) = u(x)u(y)u(z), 

where u(x) = | — x 2 for \x\ <\, u(x) = ^(| — \x\) 2 for \ < \x\ < |, and 
elsewhere (the same kernel is used in the original AP 3 M). 
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CPPA, Q = 64 



AP 3 M, Q = 64 




20 40 60 

CPPA, Q = 32 





Fig. 1. 1-D sections through initial conditions for Zel'dovich pancake test on 64 grid. 
Wavelength downwards: 64, 32, 16, 8 and 4 grid cells. Left: CPPA, right: AP 3 M. Solid line: 
density profile (arbitrary scale) dashed: velocity profile (arbitrary scale). 

A good quantitative measure of the departure of the obtained density field from 
the predicted theoretical value is the deviation of its variance. Integrating g 2 over 
x one can obtain the density dispersion given by 



\/«c - a 2 



1. 



(19) 



We obtain our results from a rectangular grid with cell size comparable to the re- 
garded wavelengths (formally we are using a cubic top-hat filter). However, the 
predicted analytical value of (S 2 )? , given by Eq. (19), ought to be corrected for 
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Fig. 2. 1-D sections through density profiles for Zel'dovich pancake test on 64 grid. 

Horizontal axis: position, vertical axis: g/(g) ■ Wavelength downwards: 64, 32, 16, 8 and 4 grid 

cells. Left: CPPA, right: AP 3 M. Solid line: analytical results, points: simulation. 

the effects of finite grid-cell size. Suppose our initial perturbation wavelength is Q 
times larger than a grid cell. The value of the predicted density dispersion on the 
grid is then given by 



{s 2 



Q 



-, Q-i 
^ p=0 



— J^ Q g(x,a)dx 



- 1. 



(20) 



In the weakly nonlinear regime the difference between Eq. (19) and Eq. (20) is 
small for long wavelengths, but it becomes important for waves with length of 
several grid cells (see Fig. 3f). 
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Fig. 3. Evolution of density dispersion in the Zel'dovich pancake test on 64 3 grid. 
Horizontal axis: scale factor a, vertical axis: (S 2 ) 1 ^ 2 a-e: simulation results for wavelengths of 4, 
8, 16, 32 and 64 grid cells. Solid line: CPPA, crosses: AP 3 M, dashed line: theoretical prediction, 
f: comparison of theoretical values given by Eq. (19) (solid line) to Eq. (20) for Q = 64 Q = 32 
(short dash), Q = 16 (long dash), Q = 8 (dot - short dash), and Q = 4 (dot - long dash). 



Results of this test for different wavelengths are plotted in Figs 3a-3e. As 
one can see, the evolution of a single wave is qualitatively correct in CPPA, but 
(especially for shorter waves) somewhat slower than predicted. The main reason for 
this is that at the final stage most of the mass is contained within walls one cell thick, 
i.e. , it has reached the resolution limit. The density dispersion curves generated with 
Couchman's code follow the theoretical curve closely until the particles come very 
close one to another. Then the accelerations become very high, the time step gets 
very short and oscillations of the wall due to interpenetration of particle streams 
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may start. 

Upon comparing we can state that the CPPA code produces reasonably accurate 
solutions on scales larger than 4 grid cells. 

4.2. Scale-free Models 

Three simulations have been performed with each code. In CPPA we used 
a 128 3 grid, while in AP 3 M 128 3 particles on a 128 3 grid. The models were 
started in an fi = 1 universe from Gaussian initial conditions with scale-free spec- 
tra (power indices +1, —1 and —3), cut-off below the Nyquist frequency, as ex- 
plained in Section 2.2. For data reduction the density fields were convolved with a 
Gaussian window of FWHM = 8 cells. 

The evolution of power spectra is shown in Fig. 4. CPPA as well as AP 3 M con- 
serve the power index on scales larger than 4 grid cells even in the nonlinear regime. 

The evolution of the density field cumulants S3 = ji^yz an d £4 = - jgzjz - » 
compared to exact results obtained with the help of perturbation theory (Juszkiewicz, 
Bouchet, and Colombi 1993, Lokas et al. 1995) is presented in Fig. 5. Results from 
both codes are in reasonable agreement with each other and with theoretical predic- 
tions. The largest difference is for £4 for n = +1, however this case often causes 
discrepancies (e.g., Catelan and Moscardini (1994) in their Monte-Carlo computa- 
tions obtained a value of 9 ± 1, while the perturbation theory gives S4 = 15.9). 
Generally the agreement is better for models with lower power index, since in this 
case there is less power in small scales, on which the codes are less accurate. 

Advanced evolutionary stages of the above CPPA and AP 3 M models are shown 
in Fig. 6 in the form of 2-D slices through the density field. Additionally, we 
performed test runs with identical initial conditions for both CPPA and AP 3 M. 
The tests were made on a 64 3 grid with CPPA and with 64 3 particles on a 64 3 
grid in AP 3 M. The initial spectral indices were —3 and +1 . The runs were started 
at z = 100, with (5 2 ) 1 / 2 = 0.03. The density fields of the evolved models are 
presented in Figs 7 and 8. One can see that all large features obtained with both 
codes are very similar, while small-scale structures (important in the n = +1 runs) 
are slightly different. Also, the density field is smooth in CPPA while in AP 3 M it 
is affected by the large scatter of the mass tracers in low-density regions. 

4.3. Performance Comparison 

Timing tests have been performed for both CPPA and AP 3 M. The test runs 
were 60 steps on SUN 10/40 for 32 3 and 64 3 particles/cells, with low density 
contrasts (no mesh refinements in AP 3 M). The results are given in Table 1. As we 
can see, CPPA performs one step 2-3 times faster. It has to be mentioned that in the 
non-linear regime in dense regions there are many fast moving particles in AP 3 M, 
while in CPPA the velocity in a given cell is averaged over all mass contained in 
that cell, thus getting smaller. Hence, due to the CFL condition, a step in CPPA is 
larger than in Couchman's code, making CPPA still faster. 




Fig. 4. Evolution of power spectra on 128 3 grid. Horizontal axis: log k , vertical axis: logP(fe). 
Left: CPPA, right: AP 3 M. Solid line: initial conditions, short dash: weakly nonlinear stage 
({5 2 ) 1/2 = 0.3), long dash: nonlinear stage ( (S 2 ) 1/2 = 2). 



5. Summary 



CPPA has proven to be a fast and useful code that well reproduces the statistics 
of cosmic fields in large scales. In the examples presented here it lacks resolution 
on small scales, but presently it is running on large grids (up to 192 3 ) on SGI 
and CRAY machines. We are also working on improving its accuracy with AMR 
(Adaptive Mesh Refinement) techniques. 



Acknowledgements. This work was supported by the KBN grant 2P-304-017- 
07. The simulations are partly performed on Cray computers at the Interdisciplinary 
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Fig. 5. Evolution of skewness (S3) and kurtosis (S<t) of the density field versus density dispersion. 
128 3 grid, Gaussian filter. Horizontal axis: (S 2 ) 1 ^ 2 of smoothed density field, vertical axis: S3 
(left) and S 4 (right). 

Solid line: PPA, crosses: AP 3 M, dashed: perturbation theory. 

Fig. 6. Cosmic density field for {S 2 }? = 2. Simulation with CPPA code (left) and with AP 3 M 
(right). Initial power indices: n = — 3 (upper), n = —1 (middle), n = +1 (lower). 

Center for Mathematical and Computational Modeling in Warsaw. 



Fig. 7. Cosmic density field for n = -3 at z = 2 . Simulation with AP 3 M (left) and CPPA (right). 
Fig. 8. Cosmic density field for n = 1 at z = 2 . Simulation with AP 3 M (left) and CPPA (right). 



Vol. 46 



11 



Table 1 
Execution times. 



problem execution time [s] 

init evol total per step 



32 3 


AP 3 M 


40 


1129 


1169 


18.81 


32 3 


CPPA 


29 


371 


399 


6.18 


64 3 


AP 3 M 


261 


8822 


9083 


147.04 


64 3 


CPPA 


249 


3987 


4235 


66.44 
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